NUMERICAL MODELING OF ELASTIC WAVES 
ACROSS IMPERFECT CONTACTS. 
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Abstract. A numerical method is described for studying how elastic waves interact with im- 
perfect contacts such as fractures or glue layers existing between elastic solids. These contacts have 
been classicaly modeled by interfaces, using a simple rheological model consisting of a combination 
of normal and tangential linear springs and masses. The jump conditions satisfied by the elastic 
fields along the interfaces are called the "spring-mass conditions". By tuning the stiffness and mass 
values, it is possible to model various degrees of contact, from perfect bonding to stress-free surfaces. 
The conservation laws satisfied outside the interfaces are integrated using classical finite-difference 
schemes. The key problem arising here is how to discretize the spring-mass conditions, and how to 
insert them into a finite-difference scheme: this was the aim of the present paper. For this purpose, 
we adapted an interface method previously developed for use with perfect contacts [J. Comput. Phys. 
195 (2004) 90-116]. This numerical method also describes closely the geometry of arbitrarily-shaped 
interfaces on a uniform Cartesian grid, at negligible extra computational cost. Comparisons with 
original analytical solutions show the efficiency of this approach. 
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1. Introduction. Here it is proposed to study the propagation of mechanical 
waves in an elastic medium divided into several subdomains. The wavelengths are 
assumed to be much larger than the thickness of the contact zones between subdo- 
mains, or interphases |18j . Each interphase is replaced by a zero-thickness interface, 
where elastic fields satisfy jump conditions. In elastodynamics, the contacts between 
elastic media are usually assumed to be perfect pQ. They can therefore be modeled 
by perfect jump conditions, such as perfectly bonded, perfectly lubricated, or stress- 
free conditions. For example, perfectly bonded conditions will mean that both elastic 
displacements and normal elastic stresses are continuous across the interface at each 
time step. 

In practice, contacts are often imperfect because of the presence of microcracks 
or interstitial media in the interphase. Take, for example, fractures in the earth, 
which are filled with air or liquid, and where jumps occur in the elastic displacements 
and elastic stresses. The simplest imperfect conditions are the spring-mass condi- 
tions (which are sometimes called "linear slip displacements"): these conditions are 
realistic in the case of incident waves with very small amplitudes 17 . The spring- 
mass conditions have been extensively studied, both theoretically and experimentally 
^120112]. This approach has been applied in various disciplines, such as nonde- 
structive evaluation of materials |22J and geophysics JJ] . 

However, very few studies have dealt so far with the numerical simulation of wave 
propagation across imperfect contacts described by spring-mass conditions. To our 
knowledge, only three approaches have been proposed for this purpose. First, Gu 
et al. developed a boundary integral method which can be applied to arbitrarily- 
shaped interfaces 0; but this method requires knowing the Green's functions on 
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both sides of the interface, which complicates the study of realistic heterogeneous 
media. Secondly, a finite-element method has been proposed by Haney and Sneider 
[5]. The jump conditions are incorporated automatically here into the numerical 
scheme, via the variational formulation, and attempts have been made to perform 
numerical analysis. The main drawback of this approach is that it requires adapting 
the mesh to the interface, which increases the computational effort. Thirdly, other 
authors have approached this problem by implicitly accounting for the boundary 
conditions by using an equivalent medium, and then deriving new finite-difference 
formulas |3] ■ This approach can be applied to arbitrarily-shaped interfaces on a 
uniform Cartesian grid, but the accuracy is low, and this method involves explicitly 
changing the numerical scheme near the interface. Note that the inertial effects were 
not investigated in any of these three cases, although they may be important factors 

na- 

The aim of this paper is to describe a procedure for incorporating the spring- 
mass conditions into existing finite-difference schemes, on a regular Cartesian grid. 
The geometry of arbitrarily-shaped interfaces is properly taken into account, reduc- 
ing the unwanted diffraction classically induced by the Cartesian grid. Lastly, the 
extra computational cost is low. For this purpose, we adapted the explicit simplified 
interface method (ESIM) previously developed for dealing with perfect contacts in 
ID ^2] and 2D \T2\ . A study has also dealt with imperfect contacts in ID In 
the present study, this approach is extended to 2-D configurations. The focus here is 
on the description of imperfect contacts; to avoid additional complications, we take 
media with simple constitutive laws (elastic and isotropic media), but the procedure 
should be suitable for dealing with more realistic media. The inertial effects are taken 
into account. 

This paper is organized as follows. In section 2, the problem is stated in terms of 
the configuration, the conservation laws and the spring-mass conditions. In section 3, 
a numerical strategy is described for integrating the conservation laws on the whole 
computational domain: the same scheme is used throughout the domain, but near 
an interface, modified values of the solution are used, which implicitly account for 
the spring-mass conditions. Section 4 is the core part of this paper: it describes in 
detail how to compute the modified values. Numerical experiments arc described in 
section 5, and comparisons with original analytical solutions show the efficiency of 
the method. Note that although no rigorous mathematical proof of the validity of 
the algorithms was obtained, the results of the numerical experiments performed were 
extremely satisfactory. 

2. Problem statement. 

2.1. Configuration. Let us consider two isotropic elastic media flo and fli sep- 
arated by a stationary interface T ffigure l2~TJl . We study a two-dimensional configura- 
tion with plane strains, and adopt Cartesian coordinates x and y pointing rightward 
and upward, respectively. The interface is described by a parametric description 
(x(t), y(r)). The unit tangential vector t and the unit normal vector n are 
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where x = 4 s and y = §7- T is assumed to be sufficiently smooth to ensure that 
x(t), y(r) and their successive spatial derivatives are continuous all along F, up to a 
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given order of derivation. 

The physical parameters are the density p, the elastic speed of the compressional 
P- waves c p , and the elastic speed of the shear SV- waves c s . For the sake of simpli- 
fication, these parameters are taken to be piecewise constant; however, they may be 
discontinuous across T 



(P, C p , C s ) 



(P0> CpO, C s o) 
(pi, Cpi, C s i) 



if 
if 



(x, y) G Q Q , 
(x, y) G fii. 



(2.2) 



The elastic fields are the two components of the elastic velocity v(v\, V2) and the three 
independent components of the elastic stress tensor cr(au, u\i, 022)- The projections, 
normal and tangential to the interface, of the elastic displacement u(u\,U2), those of 
the velocity v, and those of the normal stress cr.n are denoted by 



un = u.n, 
ut = u.t, 



vn = v.n, 
vt — v.t, 



ctn = (cr.n).n, 
<tt = (<J.n).t. 



(2.3) 



Let P be a point on V ffigure l2~T)l . and t be the time. Given a function f(x, y, t), the 
limit values of / at P on both sides of T are written 



MP, t) = lim f(M, t), 

where I = 0, 1. The jump of / across T, from fio to fij., is denoted by 

[f(P, t)} = h{P, t) - f Q (P, t). 



(2.4) 



(2.5) 



2.2. Conservation laws. To study the propagation of small perturbations in 
fi; (i = 0, 1), we use a velocity-stress formulation of elastodynamic equations. Setting 



U = T (vi,V 2 , (Til, (712,(722), (2.6) 

the linearization of mechanics equations gives a first-order linear hyperbolic system 
in each subdomain 

r\ r\ r\ 

— U + A l ^-U + B l — U = 0, 1 = 0,1, (2-7) 
at ox oy 
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which is satisfied outside T. The 5x5 piecewise-constant matrices A\ and J5 ; [l = 0, 1) 
are 
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(2.8) 



2.3. Spring-mass conditions. An incident wave at the interface generates four 
other waves: a reflected P-wave, a reflected SV-wave, a transmitted P-wave, and a 
transmitted SV-wave. To pose the problem suitably, it is necessary to define four 
independent jump conditions satisfied by the fields along T ( figure l2~T|) . Perfectly 
bonded conditions are generally used for this purpose, namely 



[u N {P, t)} = 0, 

MP, t)] = o, 



[a N (P, t)} = 0, 

MP, *)] = o, 



(2.9) 



corresponding to a perfectly bonded contact between the two solids in question. 

To describe an imperfect contact, one can generalize (I2.9|l into spring-mass con- 
ditions. With the notations defined in l|2.4|l . the spring-mass conditions are 



[u N (P, t)} = —a N0 (P, t), 

Kn 



[u T (P,t)} = —a T0 (P,t), 



d 2 

[a N (P, t)} = M N ^ u N o(P, t), 



d 2 

[<r T (P, t)} =M t -^uto(P, t), 



(2.10) 



where Kn > 0, Kt > 0, Mn > 0, My > 0, are called the normal stiffness, the tangen- 
tial stiffness, the normal mass, and the tangential mass of the interface, respectively. 
The conditions (|2.1U|) are called "spring-mass conditions" because of analogies with 
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Fig. 2.2. Spring-mass Theological model of the contact. 



the equations governing the dynamics of a spring-mass system f figure l2~2)) . One basic 
underlying assumption made here is that the elastic stresses do not affect the nature 
of the contact, and hence, that Km, Kt, Mm, and Mt do not depend on t or on the 
fields. They can vary with space, and hence with the parameter r. 

The spring-mass conditions provide an easy way of describing a wide range of 
contacts between solids, from perfect contact to disconnected media. For Km — » +oo, 
Kt —* +oo, Mm = 0, and Mt = 0, we obtain the perfectly bonded conditions (|2.9|l . 
For K N -> +oo, Kt -> 0, M N = 0, and M T = 0, we obtain a T (P, t) -> 0, which 
amounts to a perfect slip with no friction. Lastly, for Km — » 0, Kt — ► 0, Mjv = 0, 
and Mt = 0, we obtain <jm{P, t) — * and <Jt(P, t) — > 0, hence cr.n(P, t) — ► 0: the 
media fio and f2i tend to have stress-free boundaries, which means that no waves are 
transmitted from one medium to the other. 

The spring-mass conditions entail an important property: the plane waves re- 
flected and transmitted by a plane interface with conditions l|2.10l) are frequency- 
dependent These waves therefore show a distorted profile that is quite different 
from the profile of the incident wave, even below the critical angle. In addition, 
even if the incident wave is spatially bounded in the direction of propagation, the 
reflected and transmitted waves are not spatially bounded: a "coda" follows each of 
these waves. Phenomena of this kind, which do not occur with perfect conditions, are 
observed experimentally (see e.g. |17| K 

By choosing appropriate values of Kn, Kt, Af/v, and Mt, it is possible to model 
realistic configurations. The spring-mass conditions 12.1011 can also be obtained quite 
rigorously in some cases; the values of Kn, Kt, Mm, and Mt will therefore depend on 
the physical and geometrical properties of the interphase. Take a plane elastic layer 
sandwiched between two homogeneous isotropic half-spaces. If the thickness of the 
intermediate layer is much smaller than the wavelength, the spring-mass conditions 
can be deduced from an asymptotic analysis of the wave propagation behavior inside 



The spring-mass conditions model (|2.1Q(I has some limitations. First, if Km < +oo 
and Mm ^ 0, or if Kt < +oo and Mt ^ 0, these conditions are asymmetrical. In 
Appendix A, we establish that the influence of this asymmetry is either null (in 
the case of reflected waves) or negligible (in that of transmitted waves) in the one- 
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dimensional context. We have checked numerically that this is also the case in the 
2-D context. Setting up symmetrical conditions would considerably complicate the 
jump conditions, in return for absolutely negligible effects. 

The second drawback follows from the first equation given by (|2.10|) : there is no 
reason why a negative jump of un should not occur, whith an absolute value greater 
than the real thickness of the interphase. Since a penetration of both faces on the 
interphase is not physically realistic, the conditions H2.1()(l are valid only in the case of 
very small perturbations. With larger perturbations, finer modeling procedures are 
required, based on nonlinear contact laws |23|. Jump conditions of this kind, which 
are currently under study, require more complex numerical methods. 

3. Time-marching. 

3.1. Numerical scheme. To integrate the hyperbolic system (|2.7|) . we intro- 
duce a uniform lattice of grid points: (rr.j, yj, t n ) — (i Ai, j Ay,n At), where Ax — 
A y are the spatial mesh sizes, and A t is the time step. The approximation U™j 
of U{xi,yj,t n ) is computed using explicit two-step, spatially-centred finite-difference 
schemes. The time-stepping of these schemes is written symbolically 

ffr 1 = (l j) es) • (3.1) 

Hi.j is a discrete operator, and E is the stencil of the scheme. The subscripts in Hij 
refer to the physical parameters at (xi, yj). See [5] for a review of the huge body of 
literature on numerical methods for conservation laws. 

The interface T is immersed in the regular meshing, so that one can distinguish 
between two sets of grid points: the regular points, where the stencil of the scheme 
involves a single medium, e.g. f^o or f2i, and the irregular points, where the stencil of 
the scheme crosses T. The distribution of irregular points along T obviously depends 
on the geometry of T and on the stencil of the scheme. At regular points, the scheme 
(|3.1I) is applied classically, as in homogeneous media. At irregular points, however, the 
scheme is modified to take the spring-mass conditions H2.10J) into account. This 
modification is carried out using an interface method, and the main aim of the present 
article is to describe this procedure, which will be presented in detail in subsequent 
sections. 

In the numerical experiments described in section 5, we use a second-order scheme: 
the Wave Propagation Algorithm (WPALG), originally developed by LeVeque in the 
field of computational fluid dynamics [TJ31- Its stencil is (i = —2. ..2, j = —2. ..2) 
and (i, j) ^ (±2, ±2). WPALG is a useful tool for dealing with linear elastic wave 
propagation, for at least three reasons. First, it involves the use of nonlinear flux lim- 
iters that prevent numerical dispersion. Secondly, this scheme reduces the numerical 
anisotropy introduced by the Cartesian grid. Thirdly, WPALG is stable in 2D up to 
CFL=1. For the convergence measurements performed in section 5 (test 1), we also 
used the standard second-order Lax-Wendroff scheme. 

Note that other schemes can be used: in particular, we have successfully combined 
the interface method with staggered schemes, such as [T^J. It should therefore be 
possible to adapt most solvers for use with the interface method described in the 
forthcoming discussion. 

3.2. Interface method. From now on, we will focus on the time-stepping pro- 
cedure near the interface. To take into account the spring-mass conditions satisfied 
along r, the scheme ()3.1|l is also applied at irregular points, but some of the numer- 
ical values used for the time-stepping procedure are changed. The use of so-called 
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Value 




Number of 


n v 


5(* + l)(*H 


-2)/2 


components of C/f 


n c 


2(* + l)(fcH 


-2) 


jump conditions 




fc(fc- l)/2 




compatibility conditions 


n q 


variable 




grid points to estimate U k 



Table 3.1 

n with indices used throughout the text (I = 0, I). 



modified values deduced from the jump conditions is the key feature of the interface 
method developed in our previous studies: the "explicit simplified interface method" 
(ESIM) ^2E3E3. At time t n , the general method used to compute modified values 
is as follows. On each side of T, one defines a smooth extension U*(x, y, t n ) of the 
exact solution on the other side. The extension U* is built satisfying the same jump 
conditions as the exact solution U . At any irregular point, the modified value is a 
numerical estimate of U* at this point. 

Let us introduce some notations. Take an irregular point M with coordinates 
(xi, yj), belonging to fii (the following discussion can easily be adapted to the case 
where M(xi, y,j) £ Oo). Let P(xp, yp) be a point on T near M, for example the 
closest orthogonal projection of M onto T (figure |3~T1) . The vector containing the 
limit values of the exact solution U(x, y, i„) and those of its spatial derivatives at P 
up to the fc-th order is denoted by 

U k = lim T ( T U d ° T U ^— T lA (3 2) 

1 m-,p,m^i V dx^Pdyl 3 '"''<9y fe /' 

where / = or 1, a = 0, k and j3 = 0, a. This vector has n v = 5(k + l)(fe + 2)/2 
components. Throughout the text, many n with indices have been used; to avoid any 
confusion, they are summed up in table 1X11 To obtain concise expressions for the fc-th 
order Taylor expansions at P of quantities at [x-i, yj), we define the 5 x n v matrix 

< = ^77-^)1 (*< - - {J ^T 1 **) . (3-3) 

where I§ is the 5x5 identity matrix, a = 0, k and (3 = 0, a. The modified value 
U} j is then defined as a numerical estimate of the smooth extension 

U*(x I: y J: t n ) = U k IJ U h . (3.4) 

Note that Uq is the limit value of the solution and its spatial derivatives on the other 
side of r with respect to M. 
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Fig. 3.1. Zoom on an irregular point M; orthogonal projection P of M onto Y. 



The modified values at all the irregular points surrounding V are computed in a 
similar way at t n . The time-stepping at an irregular point is written symbolically 

Uft 1 = H id (u^ j+3 , (i, j)6s), (3.5) 

where 

{ u i+ij+3 ^ ( x hz> yj+~^ eni > 

= (3-6) 

The time-stepping procedure is completed over the whole computational domain by 
applying (|3.1() at the regular points. It only remains now to calculate J7*'s, since Uq 
in i|3.4fl is unknown. 

4. Calculating modified values. 

4.1. Differentiation of the spring-mass conditions. In the first step towards 
calculating the modified value Q3.4p . we look for the jump conditions satisfied by t/f 
(|3.2|) . for any k. These conditions are deduced from the spring-mass conditions l|2.10|l 
satisfied by um,t and un.t- Before describing the procedure, let us introduce a new 
notation. The vector containing the limit values of the (fc + l)-th spatial derivatives 
of U(x, y, t) at P is denoted by 

— fe+l T ( d k+1 rr d k+1 r d k+1 rr 

Ui = lim ( — - T U,..., rn U,..., r - T T l 

1 M^PMen t \dx k+1 dx k+1 ~ a dy a dy k+1 

where I = or 1, a = 0, k + 1. This vector has 5(fc + 2) components. Once again, 
the point P considered and the instant t are omitted. 

For k = 0, the two equations in (|2.1Q(I that deal with the jump in the elastic 
displacement are differentiated in terms of t. The geometry of V and stiffness and 
mass values do not depend on t; since v = ^j, we obtain 

1 d d 
[v N (P, t)) = — — <j N0 (P : t), [a N {P, t)} =M N — v N0 (P : t), 
Kn ot at 

(4.2) 

1 d d 
[vt{P, t)} = — — cjto(Pi t), [a T (P, t)} = M T — v T0 (P, t). 
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The time derivatives in i|4.2[l are replaced by spatial derivatives thanks to the conser- 
vation laws l|2.7[l . We sum up the relations thus obtained in matrix terms 

Cf{ul = C° U Q + Elul (4.3) 

C° (/ = 0, 1) are 4x5 matrices; are the vectors (I3.2[) for k = (i.e. the limit values 
of (|2.6[) h- is a 4 x 10 matrix; lastly, ill is the vector (|4. 1|) for I — and k = 0. 
Matrices describe the perfectly bonded conditions 1)2. 9fl . Matrix E® describes the 
correction induced by the springs and masses in (|2.10l) . Both matrices Cf and 
depend on r, but they are independent of t; they are dealt with in greater detail in 
Appendix B. 

To compute the conditions satisfied up to k = 1, we differentiate (|4.3|) in terms 
of t and r. First, the differentiation of i|4.3[l in terms of t yields 

The time derivatives in (|4.4|l are replaced by spatial derivatives, using the conservation 
laws (|2.7|) : with the notations (|3.2() and (|4.1I) . one readily obtains (/ = 0, 1) 

-b,)ui _» o )< 

which are injected into l|4.4fl . Secondly, the differentiation of (|4.3|l in terms of r yields 



CDt/g+cgAt/o 



L/ +^_C7 . 

(4.5) 



bmce CJ? (/ = 0, 1) and £/ depend on x(t) and y(r), the chain-rule gives 



0,1, 



d_ 

dr 



d_ 

dy 



and 



0T °~ 



x'l 5 y'l 5 ) E/j, 



a; J 5 yl 5 \ ^,2 
x'/ 5 y'/ 5 



Due to the normalization of vectors n and t in i|2.1[l and (|2.3|l . special care must be 
taken with the differentiation procedure -j^ Eq in l|4.5|l (see Appendix B). From the 
previous discussion, one builds three 12 x 15 matrices Cj, C\ and D\, and one 12 x 15 



matrix E , so that 



C\U\ = (Cl + Dl) Ul + Elul 



(4.6) 



Matrices C\ describe the influence of perfectly bonded conditions. Matrices Dq and 
Eq describe the changes introduced by the springs and masses. 

By iterating a similar procedure (k — 1) times, one can find matrices such that 



C^ = (c fc + J DS) U k +E k U k +1 , 



(4.7) 
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where Cq, C\ and Dq are n c x n v matrices, with n c = 2(k + l)(fe + 2) and n v = 
5(k + l)(k + 2)/2; Eq is a n c x 5(fc + 2) matrix (in practice and as shown in section 
14.41 the last matrix is never computed). This is a tedious procedure, however, even 
with low values of k (e.g. k = 1). This task can be carried out automatically by 
developing appropriate formal calculus tools; the simulations shown in section 5 were 
obtained in this way. 



4.2. Compatibility conditions. Some components of the spatial derivatives of 
U are not independent. We set 



A + 2 M c 2 -A 2c 2 -c 2 

^ 4(A + M ) 4(4-<*)' ^ 4(A + „) 4(t$-c2)' 1 



where A and /x are the Lame coefficients. Then, the standard plane elasticity compat- 
ibility conditions of Saint- Venant 113) are, in terms of stresses, 



d 2 an d 2 o- 22 d 2 ai2 , a 2 ctu a 2 cr 22 „ , . n . 

a 2 a o + a-i -^—o q — 5 h ai » o + a 2 Q 2 = 0. (4.9) 

or Si axay ay ay 



Equation (|4.9(l is differentiated {k — 2)-times in terms of x and y, giving the n r 
k(k — l)/2 relations 



d k (T11 9 fe cr 22 d k (T12 



a x fe --? a 9 x fc -j a ^' a ^-j'- 1 a y j+1 

(4.10) 

d^ll ^ 0-22 „ . „ , „ 

Ofl ; ;— = r— 77 + a-2 i :— s r— ^ =0, K > 2, 7 = 0, K — 2. 

aa; fe ^- 2 a^ +2 aa; fe -^ 2 a^+ 2 - 



Conditions (|4.10|) are satisfied at each point in Hi (I = 0,1), especially at P (see 
figure |2~T|) . One can therefore express the limit values U k in terms of a vector V* 
with — n m independent components 



C/f = Gfvf, J = 0,1 (4.11) 



where Gf is a n v x (n„ — n m ) matrix deduced from (|4.10|l . The relation (|4.11|l is 
useful to reduce the number of components in l|4.7|l . as seen in the next subsection. 
An algorithm has been proposed in |12) for computing the non-null components of 
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G l but there was a mistake for k > 3. The correct algorithm is (1 = 0, 1) 

a = 0, (3 = 0, 
for 7 = 0, k, for 6 = 0,.. .,7 
if <5 = then for e = 1, 5 

a = a + l, 13 = 13 + 1, Gf[a,f3] = l 
if 7 + 1 and S ^ and 7 7^ <5 then 



if 7 = 2 then 


v = 0, 77 


= u, 






else if (5 


= 1 then v = 


!? 7=1, 






else if 6 


= 7 


— 1 then 


v = 1, 77 


= 0, 




else ^ = 


1, 77 


= 1, 








a = a + 


1, 


/3 = /3 + 


1, 


G\\a,(5\ 


= 1 


a = a + 


1, 


= /3 + 


1, 


G k [a,(3] 


= 1 


a = a + 


1, 


(3 = (3 + 


1, 


G k [a,f3] 


= 1 


a = a + 


1, 


(3 = (3- 


5 + i/, 


G?[a,/3] 


= a 2 






= 13 + 


2-^, 


Gf[a,/3] 


= ai 






f3 = p + 


7, 


Gf[a,/3] 


= ai 






= /3 + 


2-Tj, 


Gf[a,/3] 


= 012 


a = a + 1, 


(3 = (3- 




Gf[a,/3] 


= 1 


+ 1 and 7 = 


5 then for e = 1 , . 


.,5 




a = a + 


1, 


/3 = /3 + 


1, Gf[ 


a,P] = 1 





(4.12) 



4.3. Solution of underdetermined systems. It is now required to express 
V\ in terms of Vq. For this purpose, we inject the compatibility conditions (|4.11() 
into the jump conditions (|4.7() . Setting 



we obtain the underdetermined system 



S k 



C k + D k )G k 



S« VI = S k V k + E k U 



k+l 



(4.13) 



(4.14) 



To find the full range of solutions of (|4.14|) . a singular value decomposition (SVD) of 
S k |16) is computed. Setting (S^)^ 1 to denote the (77^ —n m ) x n c generalized inverse 
of S\, and R k Si to denote the (n v — n rn ) x (n v — n m — n c ) matrix associated with the 
kernel of S k , we obtain 



(4.15) 



where is a (n v — n rn — n c ) vector of Lagrange multipliers. A similar procedure can 
be used to express Vq in terms of V\. 
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4.4. Numerical estimates of limit fields. Consider a set B of n q grid points 
surrounding P; for practical purposes, we define B as the set of grid points whose 
distance from P is less than a given distance q. We successively examine each grid 
point in £>, writing the fc-th order Taylor expansions of U(xi,yj,t n ) at P. If [X4, yj) 
is included in f^o, we deduce from (|3.3(l and l|4.11|l that 

eBnflo, (a*, = n,t-c/S + o(Ax fc+1 ) 

= nf J G£vS + 0(Aa; fe + 1 ) 

= lit, at (i|o) J +o(A,w), 

(4.16) 

where 1 is the identity matrix (n v — n m ) x {n v — n m ), and is the null matrix 
(n v -n m ) x ). If (xi, y-j) is included in f2i, we deduce from <|3.3|) . (|4.11|) . 

and (l^7T*51t that 

fo, yi )eBnfii, U{x hyj ,t n ) = U^U'i + OiAx^ 1 ) 

= n^. g\v\ + o(Ax k+1 ) 



II*. ({s^'s^R^ 



V k 



A 



+11^ G k (s k ) 1 E k U k Q +1 + 0(Ax k+1 ). 



Relations (|4.16() and 14.17|l are summed up in matrix terms as follows 



(4.17) 



Vq \ _ / 0(Ax fc +!) 
" • W | + 7Vt^ +1 + : |. (4.18) 

A fc / \ 0(A.T fe+1 ) 

where (U n ) B refers to the set of exact values U(xi,yj,t n ) at the grid points of B, M 
is a 5 n q x (2 n v — 2 n m — n c ) matrix, and AT is a 5 n g x 5(A; + 2) matrix. The value of 
n q depends on £>; this set is chosen so that (14.18(1 is overdetermined. In view of Table 
13.11 and Table l4~Tl this means that 

n q > (2n v -2n m — n c )/5, 
> 2(fc 2 + 5fc + 3)/5. 

To ensure this inequality, 

q= (k + 0.2)Ax (4.19) 

can be used as a radius to obtain B. From now on, numerical values and exact values 
are used indiscriminately in l|4.18|l . We also remove N U and the remainder terms. 
The least-squares inverse Af" 1 of M is computed using classical techniques, such as 
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LU or SVD [T2|. The set of Lagrange multipliers A fc is not useful for computing 
modified values; the restriction M^ 1 of M _1 is therefore defined by keeping only the 
(n v — n m ) first lines of M~ l . This gives the least-squares numerical estimate 



V k = M- 1 {U n \ 



(4.20) 



4.5. Modified value. We now have all the tools required to be able to compute 
the modified value at M(xj, yj), as defined in section 3.2. The modified value at 
(x I: yj) is deduced from l|3.4|l . (|4.11|) . and l|4.20|l 



U} J = Ul J G k M- 1 (W 



(4.21) 



A similar procedure is applied at each irregular point surrounding T. All the matrices 
required to compute the modified values are recalled in Table |4~T1 with their sizes. 



lvi_cLLI lA 




( /"i m m n "t~ 


C\ 


Hjq X Tly 


jump conditions: perfect contact (|4.7(l 


Df 


TIq X Tly 


jump conditions: imperfect contact (|4.7|l 


G\ 


Tly X {jly Tl m J 


compatibility conditions l|4.11|l 


S k 


Tic X {riy Tl m ) 


(|4T3|) 


K 


(n v —n m ) x {n v —n c — n m ) 


kernel of fif 




5 x n v 


Taylor expansions lj3.3|l 


i 


{Tly Tl rn ) X {jly Tl m J 


identity lETTTfJj) 





(riy -n m ) x (n v - n c - n m ) 


null (|^TTf^l 


M 


5 X (2 Tly — 2 Tly — Tl C ) 


BTSl) 


M' 1 


{riy - n q ) x 5n q 


restriction of M _1 $~Ttfi 



Table 4.1 

Matrices used for computing the modified values (1 = 0, 1). 



4.6. Comments about the algorithm. This algorithm can easily be adapted 
to existing codes without having to change the scheme. At each time step, one only 
needs to compute a few "modified values" 14.21fl . which do not depend on the scheme 
selected. 
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Since the interface is stationary and the nature of the contacts does not vary 
with time, most of the algorithm can be built up during a preprocessing step, con- 
sisting in determining irregular points, deriving interface conditions (|4.7|) . solving 
underdetermined systems (14.15(1 . and computing the matrices involved in (|4.21|l . All 
these quantities are stored for future use. At each time step, we need only to carry 
out one matrix- vector multiplication 1(4.21(1 at each irregular point, before running 
time-stepping procedure. The matrices involved here are small: 5 x 5 n q components, 
usually with n q ps 20. In addition, the number of irregular points is much smaller than 
the number of grid points. The extra computational cost is therefore very low. The 
CPU time measurements are likely to be approximately the same as those presented 
in P] , 

The main aim of the interface method is to incorporate the jump conditions into 
finite-difference schemes. But it also gives a finer resolution of the geometries than 
the poor stair-step description induced by the Cartesian grid. The differentiations of 
spring-mass conditions involve successive derivatives of a; and y fsee section 14. 111 . In 
addition, the Taylor expansions in 1(4.16(1 . 1(4.17(1 and 14.21(1 give an information about 
the position of P inside the mesh. 

In section l3~Tl we have stated that the interface method can be adapted to stag- 
gered grid schemes. This adaptation requires the algebra in sections 3 and 4 to be 
completely rewritten (although this is quite straightforward). To give an idea of the 
changes involved, let us focus on the staggered scheme [12] > with a grid yj, t n ) for 
stresses and a grid {x i+ i, , i„+i) for velocities. Two different sets of irregular 
points now need to be stored. The jump conditions for velocities and stresses are not 
correlated (see ((4.3(1 and Appendix B); two independant systems of jump conditions 
((4.7|) are therefore written. Since the jump conditions satisfied by the velocities yield 
as many equations as unknowns, they do not require a SVD. The compatibility con- 
ditions 1(4.1 0|l are applied only to the jump conditions satisfied by the stresses, that 
constitute an underdetermined system. The modified values calculated in sections 4-4 
and 4-5 are also split in two parts: at time t n+ i, the modified values of the stresses 
are computed (based on numerical values of er at t n ) before being inserted into the 
time-marching procedure performed on the velocities. In the same way, at time t n +i, 
the modified values of the velocities are computed (based on numerical values of v at 
^n+i) before they are used to perform the time-marching procedure on the stresses. 

4.7. Open questions. Four questions relating to numerical analysis remain to 

be solved. First, what are the effects of the term NU which was neglected in 
((4.18(1 ? In ID [TT], we established that this term is of the same order as the remainder 
terms. The proof in ID was based on performing explicit calculations on MAPLE; 
this approach seems difficult to apply in the 2-D context, because it depends on the 
geometry of T near P. 

Secondly, what is the local truncation error at the irregular points obtained by 
combining the interface method with the scheme ? In the 1-D context (111 ITS] , we 
established that the following was true: if r is the order of the scheme, the local 
truncation error at irregular points is still equal to r if k > r (which means k > 2 
in the case of second-order schemes such as WPALG), when this scheme is combined 
with the interface method. Extending this result to the 2-D context would require 
answering the previous question about N U , and bounding the matrix M~ . In 
practice, the convergence measurements indicate that this property is false in the 2-D 
context (at least for r = 2): we need k = 3 to ensure second-order accuracy (see 
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section 

Thirdly, what happens in the extreme case of a homogeneous medium (that is, 
p = pi, c p o = c p i, c s0 = Cgi, K N — > +00, K T — » M N = and M T = 0) ? 

The ideal would obviously be that U* j = U^j, in order to recover the scheme in 
homogeneous medium. In 1-D configurations |11| . we have shown that it is the case, 
given simple conditions about the set B used to estimate the modified values. In a 2-D 
context, these conditions cannot be satisfied because estimating the modified values 
requires least-squares computations. However, it would be interesting to estimate the 
difference between U*j and U^j, which might lead to an optimal choice of B. 

Fourthly and lastly, the stability analysis still remains to be performed. With a 
nonlinear scheme (such as WPALG), this seems out of reach so far, but an analysis 
could perhaps be carried out with linear schemes (such as the Lax-Wendroff scheme) . 
In the absence of theoretical rules, we have numerically considered many geometrical 
configurations, physical parameters values, stiffness and mass values. With a wide 
range of parameters, no instabilities are usually observed up to the CFL limit, even 
after very long integration times (a few thousands of time steps). However, insta- 
bilities are observed in two cases. First, instabilities can increase when the physical 
parameters differ considerably between the two sides of an interface. This problem was 
previously mentioned in J2| m connection with perfect contacts. In practice, it is not 
too penalizing when dealing with realistic media: even quite large differences between 
the impedance values, such as those encountered in the case of Plexiglass-aluminium 
interfaces, yield stable computations. Secondly, the instabilities can increase at low 
values of -Kjv,t> when the two media tend to be disconnected. This may be due to the 
small denominators in 12.10fl . In this case, increasing the order k provides an efficient 
means of improving the stability limit (see subsection !5.4ll . 



5. Numerical experiments. 



5.1. Configurations. Five numerical experiments are performed here. Except 
for test 3, the physical parameters are the same on both sides of T: this enables us 
to underscore both the effects of the spring-mass conditions and the accuracy of the 
interface method, since all wave reflections and changes observed will result from the 
spring-mass conditions, described numerically by the interface method. 

Three geometrical configurations are studied: a plane interface, a circular in- 
terface, and a non-canonical object described by cubic splines. Analytical solutions 
can be obtained for the first two configurations when the spring-mass conditions are 
constant. The analytical solution of the problem with a plane wave impinging on 
a plane interface with spring-mass conditions can be quite easily calculated using 
Fourier analysis; this procedure will therefore not be described here. We have not 
found articles dealing with the more intricate case of a circular interface; since this 
analytical solution is useful to validate the algorithm, it has been described in Ap- 
pendix C. Except for the convergence measurements in Test 1, all the computations 
were performed with the WPALG. 

Apart from test 5, the computations are all initialized by a plane right-going 
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Fig. 5.1. Test 1: plane interface between identical media at t = to + 75 i Ai. 
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P-wave 



U(x,y,t Q ) = e 



Cpl 

sin 9 

Cpl 

Ai + 2/ii cos 2 (9 



2/ii sin# cos 9 



Ai + 2 m sin 2 



V 



/ho 



x cos + y sin 9 

Cpl 



(5.1) 



where e is an amplitude factor, pi = pi c 2 sl and Ai 



Pi(cf 



pi 



2 c 2 !) 



arc the Lame 



coefficients, is the angle between the direction of propagation and the x-axis, and 
to is the initial instant. Function / is a C 2 spatially-bounded sinusoid 



/(0 



sin(w c O - i sin(2w c £) 



else, 



if < £ < 



(5.2) 



where / c is the central frequency, and lo c = 2ir f c . 

We will conclude this section by making some comments on the figures, an is 
shown with a green-red palette for P- waves, and a magenta- yellow palette for SV- 
waves (in the electronic version of this article, the plates are color). The distinction 
between these waves depends on numerical estimates of div v and curl v. Most of 
the snapshots shown are accompanied by a slice; the position of a slice is denoted by 
a horizontal line on the corresponding snapshot. On each slice, the exact solution and 
the numerical solution are indicated by a solid line and by points, respectively. 

5.2. Test 1: plane interface between identical media. Here we take the 
case of a L x x L y = 400 x 400 m 2 domain, with a plane inclined interface T. The 
points (x = 168 m, y = 17 m) and (x = 241 m, y = 253 m) belong to T, and hence 
the angle between T and the x-axis is roughly equal to 72.8 degrees. The physical 
parameters are identical on both sides of T 

Pa =Pi = 1200 kg/m 3 , 
Cpo — Cpi = 2800 m/s, 

c s0 = c s i = 1400 m/s, 

which corresponds to realistic values in the case of Plexiglass. The spring-mass con- 
ditions are constant along T; the stiffness and mass values are 

K N = 10 9 kg/s 2 , K T = 10 7 kg/s 2 , 



M N = 2000 kg/m 2 , M T = 1000 kg/m 2 . 
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Scheme N x L^ error order hi error hi order 



Lax-Wendroff 


50 


6.265e-l 




4.676e3 






100 


1.760e-l 


1.832 


1.336e3 


1.807 


+ 


200 


4.755e-2 


1.888 


3.718e2 


1.845 




400 


1.228e-2 


1.953 


9.437el 


1.978 


ESIM 


800 


3.164e-3 


1.956 


2.311el 


2.030 




1600 


7.916e-4 


1.999 


5.667e0 


2.028 


WPALG 


■50 


4.570e-l 




6.259e2 






100 


1.472e-l 


1.634 


1.504e2 


2.058 


+ 


200 


4.260e-2 


1.789 


3.954el 


1.926 




400 


1.398e-2 


1.607 


1.071el 


1.884 


ESIM 


800 


4.464e-3 


1.647 


2.833e0 


1.919 




1600 


1.503e-3 


1.570 


7.369e-l 


1.943 



Table 5.1 
Convergence measurements in Test 1. 



The incident P-wave Ij5.1|l is defined by: 6 — 40 degrees, f c — 49.9 Hz, to = 0.1 s, and 
e = 2 10~ 3 . The numerical experiments are performed with N x x N y = 400 x 400 grid 
points, which amounts to 56 and 28 grid points per central wavelength for P-waves and 
SV-waves, respectively. We take CFL=0.9, and k = 2 for the extrapolations involved 
in the interface method. To initialize the computation, one also needs to compute 
reflected and transmitted P- and SV-waves, via discrete Fourier transforms (DFT). 
To obtain the required level of accuracy, we perform the DFTs on 32768 points, with 
a sampling frequency of 0.018125 Hz. At each time step, the exact solution is imposed 
at the boundaries of the computational domain. 

Figure I5TT1 shows the solution at time steps t{ = to + 75 i At (i = 0,1,2). The 
slices are performed from 

• i = 0: (x = 70 m, y = 26 m) to (x — 360 m, y — 26 m), 

• i = 1: (x = 156 m, y — 162 m) to (x = 207 m, y — 162 m), 

• i = 2: (x = 246 m, y = 198 m) to (x = 282 m, y = 198 m). 

Looking from the left to the right on the first slice (i = 0), one can successively 
observe the reflected SV-wave, the transmitted SV-wave and the transmitted P-wave. 
As mentioned in section l2~3*l these waves do not have the same sinusoidal profile as the 
incident P-wave. A "coda" occurs after the reflected and transmitted P and SV-waves: 
as in jll| . these fields are not spatially bounded in the direction of the propagation. 

The second slice (i = 1) crosses the reflected SV-wave. The agreement between 
exact and numerical values is good; one only observes a small amount of numerical 
diffusion induced by WPALG. This is also so in the case of the last slice (i = 2) across 
the transmitted SV-wave. 

In table 5-1, we show convergence measurements performed with the Lax-Wendroff 
scheme and with WPALG, by taking k = 3. This table was obtained by taking sue- 
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cessively refined meshes with a constant CFL and measuring the differences between 
the analytical and numerical values. The orders of accuracy (2 for Lax-Wendroff, 1.5 
in norm and 2 in norm Li for WPALG) are the same as in homogeneous medium. 
Similar results can be obtained by increasing k. However, for k — 2, the convergence 
orders fall below 1.4. (we recall that k — 2 sufficies to ensure second-order accuracy 
for perfect contacts in 2-D contexts J2jj and for imperfect contacts in 1-D contexts 



5.3. Test 2: circular interface between identical media. As a second ex- 
ample, we take the case of a L x x L y = 600 x 600 m 2 domain and a circular interface 
(radius a = 100 m, centred at xq = 330.5 m, yo = 300 m). The numerical experiments 
are performed on N x x N y = 600 x 600 grid points. The physical parameters are the 
same as in the previous example. The spring-mass conditions are 



Figures f5.2l and f5.3l show the solution at the initial instant to = 0.08 s, and at fcj = 
t + 100 i At (i = 1, ...,5) on a restricted domain [100,500] x [100, 500] m 2 centred on 
the middle of the computational domain. No special treatment is applied to simulate 
the wave propagation in infinite medium (such as absorbing boundary conditions or 
perfectly matched layers), but the integration times are sufficiently short to prevent 
spurious waves reflected by the edges of the computational region from appearing in 
the restricted domain (the same comment holds for the further figures). The slices 
are carried out from 



• 


i = 


(x = 


101 


m, y 


= 301 


m) to (x 


= 289 


m, y 


= 301 


m) 


• 


i = 1 


(x = 


144 


m, y 


= 300 


m) to (x 


= 243 


m, y 


= 300 


m) 


• 


i = 2 


(x = 


108 


m, y 


= 433 


m) to (x 


= 353 


m, y 


= 433 


m) 


• 


i = 3 


(x = 


339 


m, y 


= 339 


m) to (x 


= 413 


m, y 


= 339 


m) 


• 


i = 4 


(x = 


340 


m, y 


= 300 


m) to (x 


= 493 


m, y 


= 300 


m) 


• 


i = 5 


(x = 


198 


m, y 


= 264 


m) to (x 


= 454 


m, y 


= 264 


m) 



The analytical solutions are computed on Nsessei = 90 points (see (IC.16|) in Appendix 
C), 256 Fourier points, with 1-Hz sampling frequency. 

The incident P-wave is converted into transmitted and reflected P- and SV-waves 
(i = 1,2). A series of refraction-conversion events is then observed (i = 3,4,5). The 
agreement between the numerical and analytical values is excellent. 

5.4. Test 3: circular interface between different media. As a third exam- 
ple, we consider the same circular interface as in the previous test, but take it to have 
different physical parameters on both each sides of T 

( po = 2600 kg/m 3 , c p0 = 6400 m/s, c s0 = 3200 m/s, 
(p, c p , c s ) = < 

[ p x = 1200 kg/m 3 , Cpi = 2800 m/s, c sl = 1400 m/s. 

These parameters are those of aluminium (inside T) and Plexiglass (outside T). Apart 
from these physical parameters, the values used for the computations are the same as 
in section I3T51 The initial values are the same as in figure i = 0. 

Figure 1^41 gives the results of a parametric study in terms of the normal stiffness 
Km'- perfect contact (i = 0) which amounts to — > +oo, — 10 9 kg/s 2 (i = 1), 



.11.)- 




(i = 0) 
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(i = 0) 




Fig. 5.2. Test 2-a: cylindar between identical media, t = to + 100 i At. 




Fig. 5.3. Test 2-b: cylindar between identical media at t = to + 100 i At. 




Fig. 5.4. Test 3: cylindar between different media at t = to 
normal stifness. 



+ 360 A t, with various values of 
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and Kn — 10 8 kg/s 2 (i = 2). The inertial effects and the tangential stiffness are not 
taken into account here 

Kt — ► +00, 

M N = M T = kg/m 2 . 

The computations are shown in this figure after 350 time steps. The slices are per- 
formed from x — 100 m to x = 500 m, at y = 300 m. Compared with the case of 
perfect contact (i = 0), one can see in (i = 1,2) that the signs of the reflected waves 
are reversed. The waves that have entirely crossed the interface also differ greatly 
between the three cases. The case (i = 2) is much more difficult to handle than the 
cases (i — 0,1) from the computational point of view. To obtain the good agreement 
shown in figure I5T41 (i — 2), it is necessary to use k — 4 (for the numerical solution), 
and 32768 Fourier points with a sampling frequency of 0.018125 Hz (for the analytical 
solution) . 

Smaller values of Kn have also been investigated (figures not shown here) . Below 
Kn = 10 7 kg/s , numerical instabilities have been observed when k = 2 (see subsec- 
tion E3Jl, whereas the computations remain stable when k > 3. Instabilities increase 
below Kn = 10 5 kg/s if k = 3. Stable computations and excellent agreement with 
the analytical solutions are obtained down to Kn = 10 4 kg/s if k — 4. 

5.5. Test 4: cubic spline between identical media. Up to now, we have 
dealt only with canonical objects with constant curvature. However, the use of the 
interface method is not restricted to such simple geometries. In this section, we study 
a more complex configuration corresponding to a cubic spline. The other parameters 
are the same as in section IB~3l except to = 0.07 s. 

Figure l5~5l shows snapshots of the solutions at times t% = to + 150 i At (i = 
0, 1, 2). We obviously have no analytical solutions to compare with the numerical 
values. However, the computations remain stable, and the wave phenomena seem to 
be realistic. 

5.6. Test 5: plane interface with variable jump conditions. In the last 
experiment, we study the case of a plane interface with variable spring- mass condi- 
tions. The interface T is placed horizontally at y = 130 m on a L x x L y = 600 x 400 
m 2 domain. The physical parameters have the same values as in section l5.2l and they 
are identical on both sides of T. The variable values of normal and tangential stiffness 
are 

10 12 kg/s 2 , on [x = m , x = 270 m], 

K N =K T =l g{r) on [x = 270 m , x = 330 m], (5.3) 

10 7 kg/s 2 , on [x = 330 m , x = 600 m], 

where g(r) is a cubic polynomial that ensures a C 1 continuity (t is the abscissa along 
r). The stiffness values (15.311 are those of almost perfectly-bonded media on the left 
part of the interface, and almost perfectly-disconnected media on the right part of the 
interface. The inertial effects are not taken into account: Mn = Mr = kg/m 2 . 
The computations are initialized by a transient explosive P-wave centred at (x s — 




Fig. 5.5. Test J^: cubic spline between identical media, at t = to + 150 i A t. 
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300 m, y s = 200 m), with an elastic potential ^1] 
$ip(x, y, t) = ($/ P * /) (x, y, t), 



$ /P (r, t) = 




where H is the Heaviside function, * denotes the time convolution, and / is the Ricker 
wavelet 

" n Jc 

The velocities and stresses of the incident wave can be deduced from Ij5.4[l thanks 
to classical elastodynamics relations 1 . The numerical experiments are performed 
on N x x N y = 600 x 400 grid points, with f c = 70 Hz, t = 0.03 s, t c = 0.01429 s, 
CFL=0.9 and k = 3. 

Figure 151)1 shows snapshots of the numerical solution at t = to + 150 i At. The 
palettes are changed in each snapshot to maintain a constant range of colours despite 
the geometrical spreading of the waves. On the left part of the second snapshot 
(i = 1), the incident wave is totally transmitted; on the right part, the wave is almost 
totally reflected. The rapid change in the stiffness values occuring near x — 300 m 
acts as a point source for cylindrical P- and SV-waves. On the last snapshot, one 
can also observe a rightward-moving surface wave, denoted by a yellow-magenta zone. 
Since both sides are almost stress-free surfaces in that place, this surface wave is a 
Rayleigh wave pQ . 

6. Conclusion. This paper deals with the numerical simulation of 2D wave 
propagation in elastic solids separated by imperfect contacts. These contacts are 
modeled by interfaces with linear jump conditions: the spring-mass conditions fre- 
quently used in geophysics [17j and nondestructive evaluation of materials |22| . The 
spring- mass conditions are incorporated into finite-difference time-domain schemes on 
a uniform Cartesian grid, via an interface method. This is an extension of previous 
studies on perfect contacts in 1-D ^H] and 2-D contexts, and on imperfect contacts 
in 1-D context JIJ . The interface method also provides a fine description of the sub- 
cell geometry of the interfaces, at a negligible extra computational cost. Comparisons 
with original analytical solutions show the great accuracy of this approach. 

Future paths of study might consist of taking into account numerically more 
realistic models for imperfect contact situations. A first improvement consists of 
describing the energy dissipation, as done, for example, in linear viscous bonding 
models [5]. A second improvement consists of extending the present approach for 
including nonlinear contact laws, which model realistic fractures and which prevent 
from interpenetration between the materials |23| . 




Fig. 5.6. Test 5: plane interface with variable contact at t = to + 150 i At. 
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Appendix A. On the asymmetry of 1-D spring-mass conditions. The 

aim of this section is to show that, in some limit cases, the asymmetry of the jump 
conditions (|2.1(J|) has a negligible influence. 

Consider the following simplified case: that of a compressional plane wave prop- 
agating through f2o normally to the plane interface T located at x = a. Due to the 
symmetries involved in the problem, we have a 1-D configuration, with parameters 
K = Kn, M = Mat, co = c p o, and ci = c p i (a similar study could be performed 
with a shear plane wave and parameters K = Kt, M = Mt, cq = c s q, and c\ — c s \). 
The normalized harmonic components of the elastic displacement u and those of the 
elastic stress a in medium Slo (with ko = uj/co) are: 

u(x,uj) = -ik e- lkoX +ik R ± e lkoX , 

a(x,uj) = -p a uj 2 e- ikoX - p uj 2 R ± e lk ° x , 

and in medium fii (with k\ = w/ci), they are: 

u{x,uj) = -i h T± e - ikix , 

a(x,u) = -p 1 u) 2 T ± e- iklX , 

where ui is the angular frequency, and i? ± and are the reflection and transmission 
coefficients it is required to determine (the superscripts ± are explained below) . Under 
jump conditions "shifted to the left" (as in the present paper), 

[u(a, uj)) — — a{aT, uj), [<7(a, u>)] = —Mu! 2 ii(a~, uj), 
K 

a simple algebraic procedure gives 

pi ci - po Cq - i lo{ — M 

R- = \ K 



. f pp C p! Cl \ 

Poc + pi ci + iuj y — \-MJ 



Muj 2 

2 p ci 1 



(A.l) 



po Co + pi Cl + I uj y — h M j 

Under jump conditions "shifted to the right", 

[u(a, uj)} = —a(a + , uj), [a(a, uj)} = -Muj 2 u(a + , uj), 

K 

we obtain in the same way 

. fpocopici \ 
pi ci - po c - i w M) 



R + = 



. fpocopici \ 
p c + pici+iu; y — hMJ 

. . { Po cq pi ci \ e 
Po c + Pi ci + iuj y — V M J 



(A.2) 
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Comparisons between (|A.1|) and (| A. 2|> show that the reflected waves are the same 
under both types of jump conditions, and that 



Mui 2 
K 



With realistic contacts ^2 E3> the stiffness and mass satisfy K ~ pc 2 / h and M ~ 
p h, where p and c are the parameters of the intermediate medium (or interphase) and 
h is its thickness. Since uj = 2ttc/X (where A is the wavelength) and h <C A (a basic 
assumption underlying the spring-mass model), we obtain 



< 1. 



T~ -T+ 













Hence, the coefficients of transmission are very close together under both types of 
jump conditions. 

Appendix B. Zero-th order matrices of spring-mass conditions. The 

matrices C° in (|4.3|l arc (/ = 0, 1) 



/ -y x 
x y 



C = 






-2x'y 



11 'l 'l 

-x y x — y 



o \ 



'2 

x z 
xy 



Setting 



01 



Ih 



Pa 



' sQ 



Kn yjx" 2 + ; 



Po 



Kt ^Jx" 2 + y' 2 



Va: 2 + y - 

Po 



M T 



4 = -± ^x' 2 + j/' 2 , 
Po 



the non-null components of Eq in i|4.3|) are 



E° Q [l,l] 


= 01 


y' 2 c 2 p0 + x 


E° [l,2] 


= 0i 


-2x'y'c 2 


E° [l,6] 


= 01 


~2x'y'c 2 


E° Q [1J] 


= 0i 


x' 2 c 2 p0 + y 


E° [2,l] 


= 02 


-2x'v) > 


E° Q [2,2] 


= 02 


x" 2 -y'A 


E° [2,6] 


= 02 


x 2 — y' 2 ] 


E° [2,7] 


= 02 


2sV) , 



' 2 (c 2 



pO 



2 c 2 



E°[3,3] = -/W, 
Bg[3,4]=^i', 

£°[3,9] = -/W, 
Eg [3, 10] = /3 3 z\ 
Eg [4, 3] = /3 4 x', 
Eg[4,4]=/? 42/ \ 
E?g[4,9]=/3 4 a;', 
Eg [4, 10] =&!/'• 



Appendix C. Exact solution for a plane wave on a circular interface 
with spring-mass conditions. 

This analytical solution is obtained in 6 steps: 
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1. Fourier-transforming the incident wave (|5.1|l : 

2. writing the elastic potentials on the basis of circular functions; 

3. expressing the elastic helds in terms of their potentials; 

4. computing the reflection and transmission coefficients from (I2.10|l : 

5. returning to Cartesian coordinates; 

6. inverting Fourier transform of the elastic fields (not described here). 

Step 1. The Fourier transform T of a function s(t), and its inverse Fourier transform 
J 7 ^ 1 , are denoted by 

1 r + °° 

S(u) = T(s(t)) = — s(t)e^ ut dt, 

27T J-OG 

(C.l) 

r + co 

s{t)=T- 1 {s{uo))= s{uj)e lut duj, 



where u> is the angular frequency. Applying a Fourier transform IjC.lll to the following 
part of the initial wave 1)5.1(1 

/ xcos6»i +ysin6»i 
$ IP {x,y,t) =f(t 



Cpl 

gives the harmonic potential of the incident P-wave 

$ IP (x,y,co) =e <=p^ /(w). (C.2) 

The Fourier transform of the bounded sinusoid (|5.2I) is 

f(w) = A=£-( * 2 + 5 - A 2 ) (e- l H u -lV (C.3) 

Step 2. The circular interface is centred on (xo, yo)- A point M(x, y) 6 f2i has polar 
coordinates (a; = xo + r cos0, y = yo + r sin</>). We set (i = 0, 1) 

k pl = —, k s i = —, s = e- ik "^ X0 cos9l+yo sia9l \ 9 = <j)-e 1 . 

Cpl C s l 

The first-kind Bessel functions J n satisfy the classical property ;14 

e- ircose = ^£„i n (-l)» cosn0J„(r), (C.4) 

n=0 

with e n = 1 if n = 0, 2 else. From (1C.2|) . (|C3fl . and (|C.4f> . we can therefore express 
the harmonic potential of the incident P-wave as 

& I p(x,y,u ) )=ASj2 £ ni n (-l) n cos n 6 J n (k pl r). (C.5) 

n=0 

To satisfy the Sommerfeld condition, the harmonic elastic potential $_rp of reflected 
P-waves and the harmonic pseudo-potential &rs = (0, 0, ^>rs) of reflected SV-waves 
are written on the basis of second-kind Hankel functions H n . To prevent singularities 
from occuring at r — 0, the harmonic potential $tp of transmitted P-waves and the 
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harmonic pseudo-potential V&ts = (OjOj^ts) of transmitted SV-waves are written 
on the basis of first-kind Bessel functions. Hence, we write 



+00 +00 
$rp = ^R p n cos n 9 H n (k pl r), $ RS = smn9H n (k sl r), 

n=0 n=0 

(C.6) 

*tp = J2 T " cos n 9 J n {k p0 r), * TS = J2 T n sinn9J n (k sQ r), 

n=0 n=0 



where i?^, i?* , T% and T£ are the unknown coefficients of reflection and transmission, 
which still remain to be determined. Note that numerical routines |16| usually provide 
Y n and Y n (with Y n = J n or H n ). Later in the text, we will need Y n ; to determine 
this last value, we use the differential equation satisfied by the Bessel and Hankel 
functions 

Y^(x) + \y^x) + Q) 2 ) Y n (x) = 0. 



Step 3. The elastic displacement u = T (u r ,ug) is deduced from $ (for P-waves) or 
from * = (0, 0, (for SV-waves) by 



u = grad$ for P-waves, it = curl* for SV-waves. (C.7) 
In cylindrical coordinates, the grad and curl operators are 



In cylindrical coordinates, the three independent components of <r are [Tl| 



, du r fu r , 1 due 
a rr = (A + 2 n) — h A 



d r \ r r dt 
'due u e 1 du r \ 

(777 y; 777rj- 

^ = (A + 2 M ) (ldue + u_A +x du r 



r 89 r I dr 



where A and [i are the Lame coefficients. From (JC.5I1 . (|C.6|I . I|C.9|) . we easily deduce 
the harmonic fields over the whole domain. The components of the harmonic incident 
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P-wave are 



+00 



ul p = A S fcpi £ n *" (- 1)™ cos n 9 j' n (k pl r) , 
^ ip _ AS^, , n nsmn6 

n=0 

+00 

a ip = AS Y,e n i n (-l) n cosn6((\ 1 +2 t i 1 )k 2 pl J„(k pl r) 

n=0 

+Ai^i J' n (k pl r) - Ai (P) 2 J n (k pl r)j , ( C - 10 ) 

+ OO / -. , \, 

a% =AS2^J2s n i n (-l) n sinn 6 I -J n (k pl r) - J n (k pl r) 

n=0 \ T T 

*% = ASY,Zn i n {-l) n cos n 9 (Ai k 2 pl j'i (k pl r) 

n=0 

+ (A 1 + 2u 1 )%i j;(fc pl r)-(A 1 + 2u 1 )(p) 2 J n (k pl r)\. 



The components of harmonic reflected P-waves are 



u r r p = k p i ^2 R-n cos n 9 H n (k pl r), 

n=0 

+00 

Ug P — — R p n sinnO H n (k p \r), 

r n=0 

a r / r =J2 R n cosn9 U^i + '2^)k 2 pl K(k pl r) + X 1 ^-H' n (k pl r)-X 1 i/„(fc pl r)j 

<^r£> = 2 Ml X! Si 11716 * ( "T-ffnCV 7 ") _ — H n( k pl r )j . 

+ 00 / 7 

*M = E R n cosn e { A i fc pi H 'n ( fc pi r ) + ( A i + 2 ^1) — <(*pir) 

n=0 ^ r 

-(Ai +2 M i) (^) 2 H n (kp ir )J . 

(C.ll) 
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The components of harmonic reflected SV- waves are 

+00 a 
n cosntf 



E It, CUS /K7 . . 

R n H n {k sl r), 

n=0 

+ 00 

u r g s = -k sl J2 R n sin n 6 H' n (k sl r), 

+ 00 / 1 7 

a r 4 = -2 t i 1 J2K™sn6(^ H n (k sl r) - H n (k sl r) 



n=0 
+00 



Ke = -^Y,K ^n9 (k 2 sl <(fc,ir) -^H n {k sl r) + H n (k sl r)) 
^ = -2Mi ^i^ncosnflf ^<(fc al r)-^iT n (fc al r)J . 

n— ^ ' 



The components of harmonic transmitted P-waves are 

+ OC 



(C.12) 



uf = k p0 T n cos n 9 j' n {k P or), 

71 = 

+OO 

uf = — ^ TP n sin n 9 J n (k p0 r) , 

r n=0 

+ 00 

Si* = T n cos n 9 ((A + 2 Mo ) fcpo X' (MO 

n=0 

+A ^° j' n (k p0 r) A (^) 2 J„(fc p0 r)) , ( CJ3 ) 
b% = 2Mo^ T ™ n sin n6> f J n {k p or) - — ^(fcpoOJ , 

n=0 ^ r r ' 

+00 (n k 1 

a ee = X! T ™ cos n ( A ° fc p0 J n ( fc pO r ) + (A + 2 Mo) j' n {k pQ r) 

n=0 ^ 

- (A + 2/i ) (J^j J n (k p0 r) 
Lastly, the components of harmonic transmitted SV-waves are 

+ 00 

Ug s = -k s0 J2 T n sin n 9 j' n (k s0 r), 

b% = -2 Mo ]T T„ s cosnfl (\ J n (k s0 r) - *2° , (C.14) 

n=0 ^ r r ' 

0% = -Mo sin n 9 ( k 2 s0 J„ (k s0 r) - — j' n (k s0 r) + J n (k s0 r) 

n=o ^ r r / 

n=0 
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Step 4. We now compute the coefficients R p , i?*, T p and T*. For this purpose, we 
deduce from the spring-mass conditions (|2.1U|I that 



(ujp + K p + K s ) (a + ,e) 

+ + a™) (a+,9) 
(K P e + K P e +K S o) (a\0) 



= (a*? + 3ft. - M N uj 2 (u tp + u* s )) (a-, 6), 
= (a% + a% - Mt to 2 (uf + &?)) (a", 0), 



(C.15) 



for all 6. Applying (|C.15|) to the fields (|C.10|) - (|C.14p yields an infinite number of linear 
systems. For computational purpose, these systems are computed up to NBessei terms, 
giving the systems (n = 0, 1, N BesS ei) 



Q n X, 



with X r 



(R p R s T p T s ) 



(C.16) 



and 



/ J' n {k pl a) 
- J n {k p i a) 



Y n = -ASe n i n {-iy 



Ai 



(Ai + 2/ii) fc pl J„ (fcpi a) H fepi J„(fc p i a) 

-J Jn{kplCl) 

2pL X n i\ J„(fc p i a) - J n (fc p i a) 



(C.17) 
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The coefficients of the matrices Q n are 



Q„[l, 1] = k p i H n (k pl a), Q n [l, 2] = - H n (k sl a), 
Q„[l; 3] = - ( k p0 J n (k p0 a) + — ((An + 2 fi ) k 2 p(j J n (k p0 a) 



N 



H — - k p o J n (k p0 a) - Ao ( - ) J n (k p oa) 



a 



Q„[l,4] = - ( -J n (k s0 a) - 2 ^°" ( J_ J n (k s0 a) - — j'(k s0 a) 
a Kn \a a 



Q n [2, 1] = -H n (k pl a), g n [2, 2] = k sl H n (k sl a), 



Qni 2 >3] = ~ ( - Jn(k p0 a) 



2 fj,Qn f 1 



— ^n(fc P o a) — >/„(fcpO a) 



Q B [2,4] = - f fc s0 J„(fc s0 a) + 
— kpo J n (k s0 a) + 



K T \a 2 
JM)_ 

n^ 2 



(kso Jn (k 



so a 



J n (k S Q a) 
Ai fc„ 



QJ3, 1] = (Ai +2 Ml ) fc^ <(fc pl a) + -i-^i <(fc pl a) - Ai Q ff„(fc pl a) 

QJ 3 , 2] = -2\i\n H n (k sl a) - H n (k sl a)^j , 

Ai 
a 



Q n [3, 3] = - ( (A + 2 no) k p0 J n (k p0 a) + ( — - Mjv w ) fc P o J„(fc P o a) 



-An 



/ 2 , 



in 2 n . 

, hMjvw - J n (k s0 a) 

\ a a 



2 /iQnk. 



.si I 



J n (k sQ a), 



Q„[4, 1] = 2/zm ( ^H n (k pl a) - '^-H n (k sl , 
Q„ [4. 2] = -/, , ( <(fc sl a) - ^ <(fc sl a) + ff n (fc sl a) 



Q»[4,3] = - 



/ 2^ n 
V a 2 



+ MtLO - Jn(kpoa) 



Q n [4, 4] = fc'o 4 (fcao a) - (— + Mr w : 

V a 



2 /x n fcpo 



J n (k p0 a) 



i / n \ z 

k s0 J n (k s o a) + fio ( - J Jn(k s oa). 

(C.18) 

Since the values of Hankel functions can be huge (typically 10 60 ), one must be careful 
when inverting (|C.16|) . Normalisation of Q n and the use of extended arithmetic are 
required for this purpose. 

Step 5. To return to Cartesian coordinates, we use the rotation formulas [H] 



/ "1 \ 

Oil 
V CT22 J 



/ COS ( 

sine 



V 



COS( 










cos 2 



sin <p cos cos 

,2 



sin 






2 sin (j> cos ( 

2 j. „:„2 



sin 

2 sin 6 cos 






sin 2 <fi 
sin cos ( 
cos 2 6 



\ ( v r \ 

Vg 

\ CTgg J 
(C.19) 
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